Estimate wind speeds from flight data

Flight data is recorded in the body axis. By considering the velocity vector it is possible to calculate the judging axis data. The judging axis is a bit like the wind axis, but it assumes that there is no wind. We do not have enought information in the flight log to estimate wind speeds directly because there is no airspeed sensor fitted.

If we knew the environmental wind we'd be able to convert the judging axis data to the wind axis.

If we knew the flight dynamics model then we could calculate wind axis data from the body axis data.

Option 1: Perform a flight in calm conditions, assume the wind is zero and use this as a means of estimating the model derivatives.

Option 2: make some assumptions about the model and wind then perform an optimisation over higher level parameters in order to satisfy these assumptions.

Start by creating a simulated set of flight data by assuming a wind field and some model derivatives. Then develop Option 2 and use the simulated flight data to test the process.

from flightanalysis import Section, Line, Loop, Manoeuvre
import plotly.express as px
from geometry import Points, Point, Transformation, Quaternion
import numpy as np
import pandas as pd
import plotly
plotly.offline.init_notebook_mode()

Below a template is created for some simple primitive elements. Templates are always created in the judging axis.

man =  Manoeuvre("test", 1, [
    Line(50), 
    Loop(50, -1), 
    Line(50), 
    Loop(50, -0.5), 
    Line(25), 
    Line(50, 1),
    Line(25)
])
judging = man.create_template(Transformation(Point(-50, 0,20), Quaternion.from_euler((np.pi,0,0))), 30.0)


from flightplotting.plots import plotsec
plotsec(judging, nmodels=10)

If the instantanious wind vector at every time instant is known the Section can be transformed to represent the local wind axis. Below a wind field is assumed to vary based on a power law with altitude and to be constant in time.

from flightanalysis.environment.wind import wind_power_law_builder

#the wind model is a function that returns a wind vector for a given height and time instant. 
# The wind_power_law_builder constructs the model based on a power law for a given heading, wind speed and exponent.
assumed_wind = [np.radians(30), 15, 0.2]
wind_model = wind_power_law_builder(assumed_wind)  # arg = [heading, speed, exponent]

#now use the wind model to calculate the instantious wind speed at every point in the Section
local_wind_vectors = wind_model(judging.gpos.z, judging.gt)

#plotting x and y wind speeds against height shows the assumed wind field
px.scatter(
    pd.DataFrame(
        np.stack([judging.gpos.z, local_wind_vectors.x, local_wind_vectors.y]).T, 
        columns=["height", "x wind", "y wind"]
    ), 
    y="height", 
    x=["x wind", "y wind"]
).update_layout(height=400, width=400)

Now create a wind axis section by superimposing the wind vectors on the judging axis section. The plot shows a the constructed section, with lines representing the local wind field at intervals. the Blue models show the judging axis Section and the red models show the wind axis Section.

import plotly.graph_objects as go
from flightplotting.traces import tiptrace, meshes, cgtrace, vectors
from flightplotting.traces import obj

wind_axis = judging.judging_to_wind(local_wind_vectors)

def nice_plot(judge=None, wind=None, body=None, vecs=None, obj=obj):

    scale = 4
    obj = obj.scale(scale)

    fig = go.Figure()

    asec = judge if judge else wind if wind else body

    fig.add_traces(tiptrace(asec, scale*1.85))
    if vecs:
        fig.add_traces(vectors(20, asec, vecs))
    
    if judge:
        fig.add_traces(meshes(20,judge, "blue", obj))
    if wind:
        fig.add_traces(meshes(20,wind, "red", obj))
    if body:
        fig.add_traces(meshes(20,body, "green", obj))

    fig.update_layout(
                scene=dict(
                aspectmode='data',
                xaxis=dict(visible=True, showticklabels=True),
                yaxis=dict(visible=True, showticklabels=True),
                zaxis=dict(visible=True, showticklabels=True)
            ), height=800)
    return fig

nice_plot(judging, wind_axis, vecs=local_wind_vectors).show()

Now if we make some assumptions about the FD model we can convert the wind axis data to body axis data. For now we will just consider dcz/dalpha and dcy/dbeta. This means alpha can be calculated from the lift force and beta from the side force.

Cl = -wind_axis.gbacc.z * 4.5 / (0.5 * 1.225 * wind_axis.gbvel.x**2 * 0.6)
Csf = -wind_axis.gbacc.y * 4.5 / (0.5 * 1.225 * wind_axis.gbvel.x**2 * 0.6)

dclda = 0.1
dcsfdb = 0.05

alpha = np.radians((1 / dclda) * Cl)
beta = np.radians((1 / dcsfdb) * Csf)

body_axis = wind_axis.wind_to_body(alpha, beta)

nice_plot(wind=wind_axis, body=body_axis, vecs=local_wind_vectors).show()

The principal for wind esitmation is shown below. When the wind estimate is correct, cz vs alpha and cy vs beta are linear. When the wind estimation is not correct (zero wind assumed) there is a lot of error in alpha and beta. This is illustrated by plotting the curves for three wind guesses. Guess 0 is the correct wind, guesses 1 and 2 are incorrect guesses.

#first, get the judging axis data
judge_est = body_axis.to_judging()


#make some guesses at the wind model:
wind_guesses = [assumed_wind, [np.pi, 5, 0.4], [0.4, 15.0, 0.2]]


def get_coef_data(wind_guess):
    wind_model = wind_power_law_builder(wind_guess)
    wind_vectors = wind_model(judge_est.gpos.z, judge_est.gt)
    wind_axis = judging.judging_to_wind(wind_vectors)
    alpha, beta = np.degrees(body_axis.measure_aoa(wind_axis))

    coeffs = body_axis.measure_coefficients(4.5, 0.5 * 1.225 * abs(body_axis.gbvel)**2, 0.6)

    return pd.DataFrame(np.stack([alpha, beta, coeffs.z, coeffs.y]).T, columns=["alpha", "beta", "cz", "cy"])


results = [get_coef_data(wind_guess) for wind_guess in wind_guesses]

fig = go.Figure()
fig2 = go.Figure()
for i, res in enumerate(results):

    fig.add_trace(go.Scatter(x=res.alpha, y=res.cz, mode="markers", name=f"guess {i}"))
    fig2.add_trace(go.Scatter(x=res.beta, y=res.cy, mode="markers", name=f"guess {i}"))

fig.show()
fig2.show()

To tell how close the curves are to being linear we'll do a linear regression and measure the sum of squared errors.

df=results[0]
def lincheck(x,y):
    fit = np.polyfit(x, y, deg=1)
    pred = fit[1] + fit[0] * x
    return np.sum(np.abs(pred - y))


[lincheck(res.alpha, res.cz) + lincheck(res.beta, res.cy) for res in results]

[7.311878751160823, 88.23167335944584, 33.91218524172959]

So the best wind estimate should have the lowest summed error. Now try minimizing this.

from scipy.optimize import minimize

def cost_fn(wind_guess):
    res = get_coef_data(wind_guess)
    return lincheck(res.alpha, res.cz) + lincheck(res.beta, res.cy)

res = minimize(cost_fn, [0.0, 3.0, 0.2], method="Powell", bounds= [(0.0, 4 * np.pi), (1.0, 20.0), (0.1, 0.4)])

print(f"wind test: {assumed_wind}")
print(f"wind estimate: {[res.x[0]%(2 * np.pi), res.x[1], res.x[2] ]}")

wind test: [0.5235987755982988, 15, 0.2]
wind estimate: [0.5129254533221612, 15.090789471823237, 0.19629331207024106]

have a look at the shape of the data:

import pandas as pd
from itertools import product
heads = np.linspace(0.0, np.pi * 2, 40)
speeds = np.linspace(5, 20.0, 8)
ais = np.array([0.1,0.2,0.4])

vars = list(product(heads, speeds, ais))
ress = [cost_fn(var) for var in vars]

df = pd.DataFrame(
    np.concatenate([np.array(vars), np.array(ress).reshape(-1,1)], axis=1),
    columns=["heading", "speed", "a", "cost"]
)

px.scatter(df.loc[df.a==0.2], x="heading", y="cost", color="speed")